#ifndef UNIQUESPLINE_HPP
#define UNIQUESPLINE_HPP

#include <vector>
#include <iostream>
#include <stdexcept>
#include <cmath>
#include <lapacke.h>

class SpecialBSpline {
private:
    int spline_degree = 2; // Fixed as quadratic B-splines
    std::vector<double> knots;         // Knot vector
    std::vector<double> coefficients;  // Coefficients a_i

    // Recursive computation of B-spline basis function B_{i, degree}(x)
    double B(int i, int degree, double x) const {
        if (degree == 0) {
            return (knots[i] <= x && x < knots[i + 1]) ? 1.0 : 0.0;
        }
        double left = (x - knots[i]) / (knots[i + degree] - knots[i]);
        double right = (knots[i + degree + 1] - x) / (knots[i + degree + 1] - knots[i + 1]);
        return left * B(i, degree - 1, x) + right * B(i + 1, degree - 1, x);
    }


    // Solve the linear system M * a = b to compute coefficients
    void solveLinearSystem(const std::vector<double>& t, const std::vector<double>& f_values, std::vector<double>& b) {
        int n = f_values.size();
        std::vector<double> M(n * n, 0.0); // Matrix M
        std::vector<int> ipiv(n);          // Pivot array
        int status;
        /*
        std::cout <<knots[spline_degree]<< std::endl;
        std::cout <<B(0,2,-5)<< std::endl;
        std::cout <<B(0,2,knots[spline_degree])<< std::endl;
        */
        // Construct matrix M (tridiagonal)
        for (int i = 0; i < n; ++i) {
            M[i * n + i] = B(i,2,t[i]-0.5); // Main diagonal
            if (i > 0) M[i * n + i - 1] = B(i-1,2,t[i]-0.5); // Sub-diagonal
            if (i < n - 1) M[i * n + i + 1] = B(i+1,2,t[i]-0.5); // Super-diagonal
        }
        M[0] = -B(1,2,t[0]-0.5)+B(0,2,t[0]-0.5);                          // First diagonal element
        M[(n - 1) * n + (n - 1)] = B(n-1,2,t[n-1]-0.5)-B(n,2,t[n-1]-0.5);      // Last diagonal element
        
        
    
        // Solve the linear system using LAPACK
        status = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, 1, M.data(), n, ipiv.data(), b.data(), 1);
        if (status != 0) {
            throw std::runtime_error("Error solving the system using LAPACK");
        }
    }

    // Apply boundary conditions to compute a_0 and a_N
    void applyBoundaryConditions(double f1, double fN) {
        coefficients[0] = 2 * f1 - coefficients[1];
        coefficients[coefficients.size() - 1] = 2 * fN - coefficients[coefficients.size() - 2];
    }

public:
    // Constructor: computes the spline coefficients
    SpecialBSpline(const std::vector<double>& x_vals, const std::vector<double>& y_vals, double f1, double fN) {
        spline_degree = 2;
        int n = x_vals.size();
        coefficients.resize(n + 2); 
        // n+2 coefficients (including boundary coefficients)
        knots.resize(n); 
        double length = x_vals[n-1]-x_vals[0];
        for (int i=0;i<n;++i){
            knots[i]=x_vals[i];
        }
        for (int i = 0; i < 2; ++i) {
            knots.insert(knots.begin(), x_vals[0] - (i + 1)*length/(n-1)  );
            knots.push_back(x_vals[n-1] +(i + 1)*length/(n-1));
        }

        // Construct the right-hand side vector b
        std::vector<double> b(n);
        b[0] = y_vals[0] - 2*f1*B(0,2,x_vals[0]-1.5);
        for (int i = 1; i < n - 1; ++i) {
            b[i] = y_vals[i];
        }
        b[n - 1] = y_vals[n - 1]-2*fN*B(0,2,x_vals[n-1]-1.5);

        // Solve for the interior coefficients
        solveLinearSystem(x_vals,y_vals, b);
        // Assign the coefficients
        for (int i = 1; i <= n; ++i) {
            coefficients[i] = b[i - 1];
        }

        // Apply boundary conditions
        applyBoundaryConditions(f1, fN);
    }

    // Evaluate the spline at a given x
    double getSplineValue(double x) const {
        double result = 0.0;
        for (int i = 0; i < coefficients.size(); ++i) {
            result += coefficients[i] * B(i, spline_degree, x+0.5);
        }
        return result;
    }
};
#endif
